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Abstract The Sydney University Stellar Interferometer (SUSI) has two beam 
combiners, i.e. the Precision Astronomical Visible Observations (PAVO) and 
the Microarcsecond University of Sydney Companion Astrometry (MUSCA) . 
The primary beam combiner, PAVO, can be operated independently and is 
typically used to measure properties of binary stars of less than 50 milliarc- 
sec (mas) separation and the angular diameters of single stars. On the other 
(-H . hand, MUSCA was recently installed and must be used in tandem with the for- 

Q^' mer. It is dedicated for microarcsecond precision narrow-angle astrometry of 

' . close binary stars. The performance evaluation and development of the data 

^ ' reduction pipeline for the new setup was assisted by an in-house computer 

^ . simulation tool developed for this and related purposes. This paper describes 

' the framework of the simulation tool, simulations carried out to evaluate the 

performance of each beam combiner and the expected astrometric precision of 
the dual beam combiner setup, both at SUSI and possible future sites. 
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1 Introduction 

A dual beam combiner setup was recently installed in SUSI. The main role of 
the new setup is to perform high precision narrow-angle astrometry of close 
binary stars. The relative position of one star on the celestial sphere with 
respect to another in a binary system can be determined by measuring the 
separation (in optical delay) of their fringe packet^ formed by an optical long 
baseline interferometer like SUSI. The accuracy of the projected separation 
of the binary star systems obtained from this method is determined by the 
uncertainty of the optical delay measurement. A more accurate measurement 
of the optical delay of a fringe packet can be made by mea s uring the phase 
delay of the fringes instead of the group delay (|Lawson et at 120001 ). 



Now, if the measurements are to be carried out from the ground then the 
position of the pair of fringe packets must first be stabilized because their 
positions are not static but constantly changing due to atmospheric turbu- 
lence. Th i s can be achieved using a techn ique called phase-referencing (PR) 
(|Colavital . [l99i IShao and ColavitaL Il992[) . n this technique, two beam com- 



biners are required. The phase delay of one fringe packet is measured accurately 
in the presence of atmospheric turbulence with one beam combiner (usually 
called the fringe tracker) and then fed-forward into another companion beam 
combiner to stabilize the position of the sar ne or another fringe pack et. This 
technique was demonstrated with PHASES (|MutersDaugh et ai[2010[) at PTI 



( Colavita et a i Il999l) m the near infra-red wavelengths which had achieved 



an astrometric precision of 35/ias (with separation less than 1" close b inaries) 
within 70 minutes of observation time (jLane and Muterspaugh . 20041 ). 



The dual beam combiner setup in SUSI is specifically designed to do the 
same (phase-referencing observations). The main beam combiner at SUSI, 
namely PAVO is used as the fringe tracker to measure the phase delay of 
the fringes of the primary star in the binary system in real time and the 
companion beam combiner, namely MUSCA, is used to simultaneously record 
either the fringes of the primary or the secondary star. This setup is similar to 
PHASES where both beam combiners receive the same pair of starlight beams 
from the siderostats and observe the same field of view (< 2") of the sky. How- 
ever in many other ways it is different. Firstly, PAVO and MUSCA operate in 
the visible wavelengths. Secondly, each beam combiner operates at a slightly 
different bandwidth compared to the other. Thirdly, the phase-referencing of 
stellar fringes are carried out in post-processing which eliminates the need for 
a feedback servo loop in MUSCA. Lastly, MUSCA observes only one stellar 
fringe packet at a time but can switch between a pair of fringe packets of a 
binary star during observation. 

With the introduction of MUSCA into SUSI, the existing data reduction 
pipeline was also upgraded to support the dual beam combiner configuration. 
The software development, which mainly involved putting in additional fea- 



^ a fringe packet is an interference pattern produced by a light source of finite bandwidth 
(e.g. starlight) whereby the fringe visibility diminishes quickly to zero as the optical path 
difference (OPD) producing the fringes deviates away from zero. 
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simulators code-under-test 
(this work) 

Fig. 1: Logical flow diagram of the data reduction pipeline test bench. 



tures to estimate phase delay of stellar fringes and to carry out a non-real 
time phase-referencing operation, was greatly assisted by an in-house com- 
puter simulation framework. The framework and its usage are described in 
this paper. Firstly, Section [2] gives a general overview of the simulation frame- 
work. Subsequently, Sections [31 13] and [5] describe models of fringes employed 
by the simulators to generate the test data sets while Section [5] describes a 
method that was used to include the eff'ect of atmospheric turbulence in the 
simulation. Lastly Section [7] shows the output of several dual beam combiner 
simulations and the expected performance of the instruments. 



2 Simulators and the framework 

Two simulators, which are computer models of the PAVO and MUSCA beam 
combiners, are developed to generate a set of simulated interferograms of each 
beam combiner based on user-specified inputs. The simulators were written in 
the Interactive Data Language (IDL) but the design concepts and algorithms 
described here can be implemented in any other languages. Both simulators 
can read the same set of user-specified inputs and by doing so allow users to 
simulate a dual beam combiner operation in which stellar fringes are recorded 
by the actual instruments simultaneously in real time. The simulated interfer- 
ograms can then be used to test the fringe visibility squared (V^) estimation 
and the phase-re ferencing algorithm of the upgraded data reduction pipeline 



(IKok et all . 120121 ). By comparing the user-specified input and the simulator 
output, especially the estimated since it is the main science observable of 
PAVO, the accurateness of the estimation and the performance of the dual 
beam combiner setup can be assessed. Fig. [T] illustrates the data reduction 
pipeline test bench which shows the usage of the two simulators developed in 
this work. 

An overall logical flow of both simulators is illustrated in Fig. [51 Simulation 
begins with a generic model of two or more pupils. They are then customized 
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Fig. 2: Logical flow diagram of the PAVO and MUSCA simulators. 



according to the optics of individual beam combiner. Additional and optional 
phase noise to simulate the effect of atmospheric turbulence can be included 
before the customization of the pupils. Then user-specified inputs are pro- 
cessed, e.g. to determine the number of interferograms (referred to as frames 
in the flow diagram) to be generated or the visibility and phase delay of the 
fringes to be simulated. After generating the required frame either by coherent 
or incoherent combination of the pupils, depending on the type of frame (refer 
Section . e.g. science (S-), foreground (F-), ratio (R-) or dark (D-) frame, 
the amplitude of the combined pupil is converted into photon counts and sub- 
sequently into the detector read-out units (ADU). The detector read-out noise 
is also included into the simulator output. The simulated data is then saved 
into a file of appropriate format. The simulation finishes when all the required 
frames are generated. The details of each stage of the logical flow are discussed 
in the subsequent sections. 

3 The PAVO simulator 

The PAVO beam combiner is a multi-axially aligned Fizeau-type interferom- 
eter. But unlike a typical Fizeau interferometer, PAVO forms spatially mod- 
ulated interference fringes in the pupil plane of the interferometer and then 
spectrally disperses the fringes with an integral field unit. It also employs spa- 
tial filtering in its image plane and an array of cylindrical lenslet to utilize the 
full multi-ro aperture of the siderostats at SUSI. The lenslet array fragments 
the pupil of the siderostats into several segments so that fringes from different 
parts of the pupil can be measured separately. The schematic diagram of the 
PAVO beam combiner is shown in Fig. [31 It combines starlight beams from any 
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Fig. 3: Schematic diagram of the PAVO beam combiner at SUSI. 



2 of the 11 siderostats ( Davis et al 19991) at one time. Despite that, the PAVO 
simulator developed in this work is able to simulate higher order beam com- 
bination (i.e. 3 and more beams simultaneously) because the optical design of 
PAVO at SUSI is an adapt ation of a twin instrum ent at the CHARA array on 
Mount Wilson, California (jMcAlister et a which has the capability of 

combining starlight beams from up to 3 telescopes. The original optical design 
of PAVO at CHARA and the mo dified version at SUSI h ave been discussed in 
detail bv llreland et al ( 2008f) and Robertson et al ( 2010l) respectively. 



3.1 Input and output of the simulator 



There are more than 20 input parameters to the PAVO simulator with the 
more critical ones listed in Table [TJ Some input parameters are common to 
the MUSCA simulator which will be described in Section The main output 
of the simulator is a set of FIT^ formatted files which contain images of 
simulated interferograms, hereinafter referred to as frames, as well as header 
information (e.g. timestamp, fringe lock status, etc) as would be recorded by 
the PAVO camera. There are four different types of frames recorded by the 
camera during actual observations. There are the science, ratio, foreground and 
dark frames. Examples of each of these frames are shown in Fig. 2] Each frame, 
except the dark, contains images of the spectrally dispersed (horizontally in 
the figure) pupil as sampled by an array of lenslets. The number of lenslets, 
A^LL, is different between PAVO at SUSI and at CHARA. The example shown 
in Fig. m is of the former which has the lenslets arranged in a one- dimensional 
array and there are 4 lenslets per pupil. The left pupil is for science while the 
middle and right pupils are used for tip-tilt correction and therefore are ignored 
in the simulation. Frames of PAVO at CHARA contains only the science pupil. 



3.2 Model pupil 



Before the pupils are combined to form fringes, they are first spatially filtered 
in the image plane with square apertures (hereinafter referred to as the PAVO 

^ |http : //f its .gsfc .nasa.gov I 
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Tabic 1: Input parameters for the PAVO (P) and MUSCA (M) simulators 



Name 


Simulators 


Description 


iSTART 


P,M 


Start of simulation in Julian date 


tSTEP 


P,M 


Exposure time of camera/photodetector 


-^S-FITS 


P 


Number of science type FITS to generate 


-'Vr-FITS 


P 


Number of ratio type FITS to generate 


A^F-PITS 


P 


Number of foreground type FITS to generate 


-'Vd-FITS 


P 


Number of daik type t'iiS to generate 


-'Vmed 


P,M 


Number and types of optical media 


c 


P,M 


Astrometnc UPU m m 


A^TEL 


P,M 


Number of telescopes 


I > 
I > 


r ,iVi 


Details of telescopes (e.g. baselines) to be used 


my 


P,M 


Magnitude of source in V band 


V 


P,M 


Model complex visibility of source 


D,d 


P,M 


Offset of i| in m 


ro 


P,M 


Fried parameter in m 


-1 

O-ro 


P,M 


Wavelength in which ro is specified (fim) 


TO 


P,M 


Coherence time in milliseconds 


Lo 


P,M 


Outer scale of atmospheric turbulence in m 


iSTEP 


M 


Time interval between steps 


A^gxEP 


M 


Number of steps per scan 


A^SCAN 


M 


Number of scan to simulate 


-^^SCAN 


M 


Length of a scan in fira 



mask), one for each beam, to remove high spatial frequency noise arising as a 
consequence of atmospheric turbulence and optical aberration. A pupil from 
a telescope is modeled by a square matrix, P, of size A^ppx X TVpFT; 

Jnfft.Nfft ; perfect wavefront 

Jnfft,Nfft ° "^xp ; corrugated wavefront 

where J represents a unit matrix consisting of all Is of size (row x column) 
indicated by its subscript and the phase component is either zero or (p so as 
to represent a perfect or corrugated wavefront. A corrugated wavefront due to 
atmospheric turbulence will be elaborated in Section [S] The notation o is the 
Hadamard or the element- by- element multiplication operator of two matrices. 
The model of a spatially filtered pupil, P', is then, 

p' = -Fr"i{no^r{o°p}} (2) 

where ^-^d □ are square matrices of the same size which define a circular 
pupil and the spatial filter respectively. Each element in the matrices is defined 



if ^{u - 7Vfft/2)2 + iv- A^fft/2)2) < iVDiA/2 
otherwise 

if \x - Nfft/2\ < A^mask/2 and \y - Nfft/2\ < A^mask/2 
otherwise 

(3) 
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Fig. 4: (a) An example of real on-sky data of a Gru recorded by the PAVO 
camera. It contains images of spectrally dispersed (horizontally) pupils as sam- 
pled by an array of four lenslets. The science pupil on the left (indicated by 
the dotted box) of the frame is where fringes are formed by a pair of over- 
lapped pupils. The pupils in the middle and left of the frame are used for 
tip-tilt control. They are ignored and will not be generated during simulation. 
The images of the science pupil for different frame types are shown in (b) and 
(c). The images, from left to right and top to bottom, belong to four different 
types of frames, namely science, foreground, ratio and dark respectively. The 
images in (b) are recorded by the camera while those in (c) are generated by 
the simulator. The convention adopted by the mathematics in this paper takes 
the vertical axis, which is the fringes direction, as the x-axis and the horizontal 
axis, which is the spectral channel direction, as the y-axis. 



where iVoiA and A^mask are relative sizes of the diameter of the telescopes 
and the width of the PAVO mask respectively. Afft, A^dia and A^mask take 
only integer values. The ratio of the width of the mask to the size of the image 
varies with the wavelength of light and is defined in Eq. (0]) . In this paper the 
wavelength of light is always represented by its reciprocal, or wavenumber, a or 
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cTj , where j is the index of vector cr which represents a range of wavenumbers 
appUcable to PAVO. The wavelength scahng factor in Eq. (U) is derived based 
on the /-ratio of the beam and the distance of the PAVO mask from the 
pupil plane. The different values between PAVO at SUSI and CHARA means 
the optical setup at the two interferometers are not exactly the same. The 
simulator keeps the denominator of the left-hand side (LHS) of Eq. constant 
and varies A'mask at different wavelengths of light to satisfy the equation. 

A^MASK _ r 1.7 a J ;SUSI 
Nfft/Nbia 1 8.5ctj ; CHARA ^ > 

With A^FFT and A^dia typically set to 256 and 114 respectively, the range of 
values of A^mask is 5-7 for SUSI and 22-30 for CHARA. 

Before a model pupil is used to form fringes, it is down-sampled to the size 
of an actual image taken by the PAVO camera, which is A'x x A^ll, where 
Ax is the number of pixels in the x-axis and All is the number of lenslets in 
the y-axis. The values of each parameter are listed in Tabled As a result, the 
model pupil becomes, 

^ ^ AT AT -^NffT,Nx ^NffT.Nll (5) 



where Xm.n is a MxN matrix. 



'M,N — 



Jm/n,i 



J 



M/Na 



(6) 



J 



M/N,l. 



and Jm/n,i is a M/N x 1 unit matrix consisting of all Is. This simple down- 
sampling only works if M is an integer multiple of A. For example, to down- 
sample a 6x6 matrix P' to a 2x3 matrix P, the following operation can be 
applied, 





"1 0" 


J. 




1 




1 


1 




6 


1 






1 






1 





p' p' . . . p' 



P' 

2,1 



P' 



pt pi 



P' 

2,6 



P' 





"1 0" 




1 




1 




1 




1 




1 



(7) 



where the first and the last matrices on the right-hand side (RHS) are Xe.2 
and X6,3 respectively. 

The values of the physical parameters in Eq. Eq. (O and other equa- 
tions in this section are listed in Table [21 
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Table 2: Specification of PAVO 



Parameter 


Notation 


PAVOOCHARA 


PAVOO 


iSUSI 


Number of spectral channels 


N„ 


19 


21 




Spectral range (^m~^) 


a 


1/0.88 - 1/0.63 


1/0.80 - 


1/0.53 


Number of lenslets 


Nll 


16 


4 




Number of pixels in x-axis 


ATx 


128 


32 




FOV to pupil size ratio 


R 


5.8 


1.8 





3.3 Simulating various types of frames 

Different types of frames are generated using different combinations of model 
pupils described in Section 13.21 Sets of 50 frames are saved into a FITS file. 
Each FITS file has a timestamp of, 

ta = tSTART + 50 X QfisTEP (8) 

where a is the number of FITS files already generated before the current one. 
The typical value of isTEP for PAVO at SUSI and CHARA is 5ms and 8ms 
respectively. The number of files to be generated for each type of frame is 
determined by the user. 

3.3.1 Science frames 

The science frames are generated using two or three pupils, depending on the 
number of telescopes in use. Simulation of PAVO at CHARA can use up to 
three. In reality the pupils are aligned to overlap each other and combined 
to produce spatially modulated fringes across the pupils. The model intensity 
across the overlapped pupils at one particular wavelength is given as, 

F = w,J2 E ^e,e aeP^'^ ° ^e^^'^ ° exp (9) 

where the notation P represents the complex conjugate of the variable P. 
The indices and 6 denote one of the several pairs of telescopes used in the 
simulation (A^tel) while ao denotes the weighted amplitude of P^^^ such that, 

al+al + ---+ a%^^^ = 1 (10) 

This condition has no physical reason but is imposed for the convenience of 
scaling the normalize intensity to the right photon rate in Section [5] The term 
Wj states the relative intensity of the summation in Eq. ^ at one wavelength 
while the vector w which Wj is a part of describes the spectrum of the light 
source and the bandpass profile of PAVO, 



w = [wi W2 ■ ■ ■ wn„] ; where wJn„,i = 1 



(11) 
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This term can easily be customized by user according to the need of a simula- 
tion. However, the results shown in Section [7] were simulated with a smooth- 
edged top hat function for w. The model of fringes, F, is expressed in this 
form in order to allow a model of complex fringe visibility, Vg g, to be applied 
to the pairs of pupils. The complex fringe visibility matrix, which is a user 
supplied input, is defined as. 



V = 



V{a,B,,2) 



TVt 



2. N'T 



^^(o-j52,7Vtel) 
1 



(12) 



where each off diagonal element represents the fringe visibility of a model light 
source at a given wavelength and baseline. The term Eg g is the magnitude of 

a baseline vector Bg g formed by a pair of telescopes 6 and 9. Lastly the matrix 
# represents an additional phase difference between the two pupils and it is 
used to model the difference in piston and tilt in the wavefront of the pupils. 



N.iSg si?xJi, 



Nr 



(13) 



The first three terms in Eq. (|13p represent the piston term. The first term 
(g g is the astrometric OPD due to the position of the target star with respect 

to the baseline Bg g while the second term [NJj.^z^^'^^ is the optical path of 
the delay line used to compensate the astrometric OPD. 



c = 







S-Bi.2 





-s • Bi, 



B 



2,Ni 



S ■ Bi^jVtbl 
S ■ B2,7Vtel 







(14) 



The optical delay line can comprise of various optical media. The types of 
optical media are specified by the user but practically it is not more than 
4 different types (e.g. vacuum, air and two types of glass, BK7 or F7). The 
refractive indices for each medium at the w avenurnbers o f PAVO are calculated 
using values and constants obtained from iTango (|l990t ). The notation [N]j,, 
represents the j-th row of the refractive indices matrix, N, where. 



N 



ni((Ti) 
"1(0-2) 



"2(0-1) 
"-2(0-2) 



"A'med 
"-^Vmed 



('^l) 



"i(o-Ar^) "2(o-Ar^) ••■ "ATmed (o-^V,, ) 



(15) 



and rii is the refractive index of one optical medium. In order to set a con- 
vention, the first medium (i = 1) is air. Each element in the column vector z 
represents the optical path length of each medium and to optimally compen- 
sate a given astrometric OPD the values of z are calculated using the method 
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described bv iTango ( 1990h . The third term in Eq. is an offset term to 
allow users to simulate a non-optimally compensated astrometric OPD. The 
user input matrix D is defined as, 



D 





-di.2 



'di. 



di,2 




J-2.Ni 



dl.N-TEl, 







(16) 



where di^2 for example is the OPD offset for baseline Bi 2- 

The last term in Eq. (jl3p is the OPD caused by the differential wavefront 
tilt between a pair of pupils at the pupil plane. It is proportional to the sep- 
aration of the apertures on the PAVO mask and inversely proportional to the 
distance of the pupil plane from the mask. Sg g is the ratio of the separation of 
the apertures on the PAVO mask to the width of each aperture and is defined 
as. 



216'- 



(17) 



R is the ratio of the field of view (FOV) of the PAVO camera to the diameter 
of one pupil at cr = l^m^^ and its value is given in Table [21 Lastly, x € 
M : — 1/2 < < 1/2 is a column vector of length A'x which represents the 
pixels across the field of view of the camera along the direction of the tilt. 
This direction is also the axis where interference fringes are formed across the 
camera and is referred to as the a:-axis by convention. 

Now, F is defined at just one wavelength. In reality, the combined pupils are 
spectrally dispersed by a prism. In order to model this F is evaluated N^, times, 
each time at a different wavelength within the PAVO spectral bandwidth. 
Multiple F matrices are then rearranged in the following order to mimic the 
actual interferogram recorded by the camera using only their real parts (as 
denoted by the notation 5R), 



PHOTONS 



*,0 



OR 



(1) 



*,1 



OK 



(1) 



} (18) 



+ e 



The superscripts in parentheses represent the matrices evaluated at differ- 
ent wavelengths within the spectral bandwidth. The interferogram, which 
takes only the real part of the F, is padded with columns of zeros, 0. This 
is done according to a PAVO parameter definition file. The file describes the 
pixels locations of a spectral channel within the camera's field of view which 
was determined through calibration with up to two lasers. The scaling factor, 
-/VpHOTONS, in Eq. converts intensity to energy in terms of number of 
photons. This factor is proportional to the brightness of the target star, num- 
ber of telescopes used and the exposure time, all given by the user. A noise 
term, e, is added to the interferogram to simulate photon noise, multiplication 
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noise and read noise of the EMCCD camera. It is not purely an additive term 
as suggested because the expression of the noise term in Eq. (ITSl) is simphs- 
tic. Physical models are used in the simulation to generate the photon and 
multiplication noise components based on the number of photons. 

3.3.2 Foreground frames 

Foreground frames are generated in a very similar way to the science frames. 
Instead of setting D to values within the coherent length of the PAVO spectral 
bandwidth, it is set to a very large number (e.g. Im) so that no fringes are 
formed across the pupils. Furthermore the visibility of the fringes are set to 
zero. In the actual beam combiner, foreground frames are recorded by giving 
the optical delay line a large offset from its last position where fringes were 
found. 

3.3.3 Ratio frames 

Ratio frames are generated using only one pupil at a time. In the actual beam 
combiner, such frames are recorded when one of the many beams is blocked 
from reaching the camera. This type of frame is used by the data reduction 
pipeline to determine the intensity of pupil from each telescope. With only one 
pupil, the F matrix becomes, 

F = aePC') (19) 

3.3.4 Dark frames 

Unlike previous types, dark frames are generated without any pupil. In the ac- 
tual beam combiner, dark frames are recorded when all the beams are blocked 
from reaching the camera. The interferogram contains only the noise term 
which in this case made up of only the read noise of the camera. This type of 
frame is used by the data reduction pipeline to subtract the noise floor in the 
interferogram. 

J^ = e (20) 



4 The MUSCA simulator 

The MUSCA beam combiner is a co-axially aligned pupil-plane Michelson 
interferometer. It combines only two beams at one time, each from one sidero- 
stat. There are no spatial fringes in the image plane. If the image is diffraction 
limited, the Airy disk will be completely dark when the pupils are out of phase 
and bright when the pupils are in phase. MUSCA produces interference fringes 
by varying the difference in piston between the two pupils through time. This 
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Fig. 5: Schematic diagram of the MUSCA beam combiner at SUSI. The beam- 
sphtter in the diagram is labeled BS while the APDs are single pixel avalanche 
photodiodes used for photon-counting. 



Table 3: Specification of MUSCA 



Parameter 


Notation 


Typical values 


Spectral range (/im~^) 


<T 


1/1.0 - 1/0.77 


Scan range (fim) 


-t'SCAN 


30, 140 


Number of steps per scan 


A^STEP 


256, 1024 


Time interval between steps (ms) 


*STEP 


0.3 



is done by changing the optical path length of one pupil in air by moving a mir- 
ror in discrete steps between two locations back and forth rapidly. The mirror, 
hereinafter referred to as the scanning mirror (N18M in Fig. [5]), makes a scan 
by moving from one extreme position to another. The physical parameters 
related to the scanning mirror and the operational spectral range of MUSCA 
are listed in Table [H The schematic diagram of the MUSCA beam combiner is 
shown in Fig. [SJ Fringes are formed at the beamsplitter (BS) and are recorded 
by a pair of avalanche photodiodes (APDs), one on each s ides. However th e 
details of its optical design are discussed in another paper (|Kok et a 



4.1 Input and output of the simulator 

The format of the input to the MUSCA simulator is exactly the same as the 
input to the PAVO simulator. Some input parameters are common to both 
simulators but there are some parameters which are applicable only to the 
MUSCA simulator. The parameters are listed in Tabic [TJ 

Instead of FITS, however, the output of the MUSCA simulator is a plain 
text file, which has the same format as one generated by the actual beam 
combiner. The file contains a time series of photon counts recorded by the 
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OPD (|im) 

(a) 

Fig. 6: An example of two sets of time series of photon counts (a) recorded by 
the actual MUSCA beam combiner and (b) generated by the simulator plotted 
as separate scans. Each scan is plotted with an offset for visualization. The 
cc-axis shows the relative OPD within the scan range of the scanning mirror. 
The top and bottom plots are scans of different direction. 



APDs in the image plane as the scanning mirror periodically scans through a 
predetermined scan range. Each photon count has a timestamp with a precision 
of 10 microseconds. Fig.[B]shows an example of the photon counts recorded by 
the actual instrument as well as a set generated by the simulator. 



4.2 Model pupil 



The model of a pupil in MUSCA is straightforward as there is no spatial filter 
in MUSCA. Suppose the model pupil is represented again by the matrix P of 
size TVpFT X A^FFT7 then, 

P = 0°P (21) 
where P and Q ^-re defined earlier in Eq. ([T|) and Eq. respectively. 



4.3 Simulating the photon counts in a scan 

The temporal fringes in MUSCA are generated using a pair of pupils, P^^^ and 
P*-^-*. In the physical instrument the pupils are aligned and combined coaxially 
at a beamsplitter to produce two sets of output pupils. In the simulator the 
output pupils are modeled as two No- x iVsTEP matrices, 

F(^) = 1 + 2aia25R {((w • v)TJi_Nstep) o Q o exp (z*)} 
f(^) = 1 - 2aia23fi {((w • v)'^Ji,Nstep) o Q o exp (z*)} 
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The matrices and their superscripts (L) and (R) represent the intensity of the 
output pupil on the left and right side of the beamsplitter. The variables ai, 02 
and w are previously defined in Eq. (ITO)) and Eq. ((TT|) . The vector v, which is 
similar to the matrix V in Eq. ([T^ . is the model fringe visibility of the source 
for the baseline defined by telescopes of pupil P*-^^ and 

V = [V^(aiBi,2) F(a2Si,2) • • • F(aw„Bi,2)] (23) 

Each element of Q is, 

EWfft y^Wfft p(l j',fe) p(2,j,fc) 
^- f_ — M^l l^V = l ^I^.V (24) 

"'' Area of circular aperture in Q 

The real part of Eq. ([M)) . di{Qj^k}, represents the normalized intensity at 
wavenumber ctj due to the sum of pupils P*^^^ and P*^^^ at fc-th step of a scan. 
$ is the additional phase difference between two telescope pupils at various 
wavelengths, 

* = 27r(crT^ + J^^io-NZ + crTd) (25) 

and it is used to model the difference in piston between the pupils at each step 
of a scan. 

The structure of Eq. (|25p is very similar to Eq. ((T^ . The first term in 
Eq. C, is the astrometric OPD per scan of the baseline B1.2. Each element, 
Cfc, is evaluated at time ta,k, 

a = Is -61,21,=,^^ (26) 

where. 



ta,k = iSTART + ^STEP ("A^STEP + k) (27) 

and a is the number of elapsed scans before the current one and k is the 
index of a step within a scan. The typical value of ^step for MUSCA is 0.3ms. 
The second term in Eq. NZ, is the optical path of the delay line used 

to compensate the astrometric OPD. The matrix, N, previously defined in 
Eq. (jl5p . denotes the refractive indices of each optical medium in the path 
while the matrix, Z, denotes the path length of each medium at every step in 
a scan. 



Zl{ta,l) Zi{ta,2) ■ 
Z2{ta,l) Z2{ta,2) ' 


^1 (ta.ATsTEP ) 
Z2{ta,Ns-T-E-p) 


+ 


Ol,NsTEP 


.^JVMED(^a,l) 2AfMED(^Q,2) • 


^A^MED (^C(,-'VSTEP ). 




_Ol,NsTEP_ 



The values of elements in the first term of t he RHS of Eq. (051) are calculated 
using the method described bv lTangol(|l990l) . The vector I in Eq. (pS)) describes 
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the relative change in the optical path length of air at every step in a scan due 
to the motion of the scanning mirror. 

^ ^ i^SCAN ^, 

^STEP - 1 (29) 
[-iVsTEp/2 -iVsTEp/2 + -lOl-- - iVsTEp/2 - 1] 

Elements in the same column of Z have the same timestamp as an element 
with the same index in ^. The last term in Eq. d, is the user-specified 

offset at each step to simulate a non-optimally compensated astrometric OPD. 

d= [di d2 ■■■ dATsTEp] (30) 

After the pupils are combined and have formed fringes, it is assumed that 
the entire image of the pupil falls within the active area of the photodiodes 
and all photons are detected. The photodetectors in MUSCA have only one 
pixel and are unable to resolve any spatial variation in intensity at the image 
plane. Therefore only an average intensity is recorded, hence the term Q in 
Eq. (j22p . In addition to that the photodetectors in MUSCA are unable to 
resolve intensity variation across wavelengths. Therefore the number of photon 
counts recorded by the photodetectors is the sum of photons across the entire 
MUSCA operating bandwidth, 

^ = iVpHOTONS X Ji,N„F + e (31) 

Similar to Eq. (jlSp . iVpHOTONS is a scaling factor that converts the intensity 
of the output pupil to the number of photons expected from the source given 
its brightness in magnitude scale and e is a noise term included to simulate 
photon noise and the dark count noise of the detector. 



5 Converting magnitude scale to photon rate 

The scaling factor iVpHOTONS in Section [3] and Section 2] is estimated based on 
the expected throughput of the beam combiner to be simulated, the efficiency 
(Q.E.) of the APDs and a calibrated magnitude-to-flux scaling factor, _Fiy, by 



Besselll ( 19791 ). The first two factors collectively describe the efficiency of the 



instrument, rj, which is found to be -3% for PAVO and -1% for MUSCA. The 
lower efficiency in MUSCA is possibly due to the aluminium coated mirrors 
used in SUSI and a silvered beamsplitter in MUSCA. The values of Fi, at 
differen t phot ometric bands are listed in Table SI which is reproduced from 
Besselll ^97^. From i^jy, the number of photons from a my magnitude star 



collected by a telescope with an area of Atd m^ in At seconds can be estimated. 
Putting the factors together, 

iVpHOTONS 77 X 1.51 X 101° X l0-"-/2-5 X F.^^^AelAt (32) 
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Table 4: Absolute flux calibration of a Lyrae ( Bessell . 1979() 



Filter 




AX 


Flux density, 


band 


(fim) 


(fim) 


(xlO-2^ W m-2 Hz-i) 


U 


0.36 


0.076* 


1.81 


B 


0.44 


0.094* 


4.26 


V 


0.55 


0.088* 


3.64 


Rc 


0.64 


0.57-0.72 


3.08 


Ic 


0.79 


0.725-0.875 


2.55 



Ihttp : //www . astro . umd . edu/- ssin/ASTR620/mags ■ html | 



where the A{l/a) is the bandwidth of the beam combiner and AX is the 
bandwidth referred from Table H] which has an effective wavelength close to 
that of the beam combiner. As an example, in the case of the MUSCA beam 
combiner, Z\(l/cr) ~ 0.29Aim, AX = 0.15^m, X^s = 0.79, Atei = 7 x IQ-^m^ 
and At ~ 80ms (256 x 0.3ms). 



6 Simulating the effect of turbulent atmosphere 



A turbulent atmosphere introduces amplitude and phase fluctuation to an 
otherwise plane wavefront of light from a distant star. Although amplitude 
fluctuation can be simulated, both PAVO and MUSCA simulators make the 
assumption that the atmospheric phase fluctuation only gives rise to pertur- 
bation in the phase of a wavefront. This approximation by discarding the 
amplitude (scintillation) term, also known as the near-field approximation, 
holds very we l l for pupils larger than 2.5cm under typical turbulence condi- 
tion (|Roddieil . lT98l . 

The phase fluctuation in the atmosphere is simulated using a large (A^atm x 
A'atm) two-di mensional array of random phasor, which has a power spec- 
trum given by (jRoddieii . Il981t iGlindemannl . 120111) , 



= C (Lo 2 + («2 



-11/6 



(33) 



The notation (f) denotes the Fourier transform of 0, u and v are the indices of 
the array and Lq is the outer scale of the turbulent structure of the model at- 
mosphere, which is set to a very large number in the simulation. The phase fluc- 
tuation is recovered by taking the inverse Fourier transform of 4'u,v exp(ie) 
where e is the phase of the Fourier transform and it is a random variable. The 
randomness of the generated phase fluctuation, (f>, is controlled by adjusting 
the scaling factor, C. It is tweaked so that the structure function of (p, 

- </'«H-.„,.+r„ I') = 6.88 (:M±H] (34) 



'^0 



has its characteristic Fried parameter, rg, set to the value specified by the user. 
Fig. [Tl^a) shows an image of the random phases generated using this method. 
The grayscale of the images indicate the value of the phases. 
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1^3 J 




Fig. 7: (a) An array of random phases generated with the inverse Fourier 
transform method. The value of the phases is indicated by the shades of gray 
in the images, (b) Small portions of the larger array in (a) sampled at slightly 
different positions at each time step showing the progress through time. Time 
increases from left to right and top to bottom. 



On the other hand the phase fluctuation of a wavefront across a telescope, 
or (p in Eq. ((T|) and Eq. (HH, is simulated by extracting a small portion 
(~ 100x100) of the much larger 4> array. The typical size of the 4> array is 
2048x2048. Since tq is specified at a certain wavelength, cr^^^, and the phase 
fluctuation in <p is generated according to rg, the value extracted from the 
larger array must be scaled by a factor of <Tj/arg, where (Tj is the desired 
wavenumber for simulation, before applying it to the model pupils in Eq. ^ 
and Eq. (|2ip . The position of this small sampling window is displaced across 
the larg er array after e very time step of the simulation. This simulates the 
effect of lTavloir s (Il938h hypothesis of a frozen atmosphere drifting across the 
aperture of the telescope. The rate of displacement of the sampling window 
depends on the wind speed which is estimated from, 



V = 0.314 



ro 
To 



(35) 



where tq is the coherence time of the phase fluctuation and the time step of the 
simulation, tsTEP- Both parameters are specified by the user. The direction of 
the displacement is random but remains constant throughout the simulation. 
Fig. [71[b) shows several snapshots of phase variation over a small portion of a 
larger array drifting across the sampling window. 

A separate array of cp is generated for each telescope. By doing this, it is 
assumed that the phase fluctuations over individual telescopes are uncorre- 
lated. As an effect, the low- frequency phase fluctuations do not increase with 
baseline in the simulation. In practice this scenario is true if the baseline is 
longer than the outer scale of tu rbulence, Lq, which is in the order of 100m in 
the troposphere (jRoddieii . 119811 ). Therefore this approach of having a separate 
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array of </> for each telescope does not simulate the piston term of an aber- 
rated pupil due to atmospheric turbulence realistically. In order to address 
this shortcoming in the PAVO and MUSCA simulator, the input variable, D 
or d (refer Table [T]) , can be used to add an extra differential piston between 
pupils from two apertures. 

Another shortcoming of this method in simulating the phase variation 
across a turbulent atmosphere is a structure function that deviates from its 
theoretical value at large distances. Fig. |5] shows comparison between a theo- 
retical structure function and one calculated from the (p. The plots in the figure 
are versus r and normalized to the size of the array, A'atm- The 

theoretical Kolmogorov (KM) model curve is a straight line with unity slope 
(rg = 1). The deviation from the theoretical value is very pronounced espe- 
cially at large distances. This is not surprising because the power spectrum, 
|(/>P, is undersampled at low frequencies, where most of the energy resides. 
A practical solution to this problem, which is the approach implemented in 
this work, is to sub-sample the phase fluctuation fr om a much larger array 
(jMcGlamervl . Il976l: IShaklanl . 1 19891 : iLane et all Il992l) . With a telescope aper- 
ture model of ~5% of the size of the larger phase array, such approximation 
produces phase fluctuation that is at most 20% off the theoretical KM model. 
This is shown in Fig. |S1 Also shown in the figure is a theoretical von Karman 
(VK) model which has a structure function very similar to that of the gener- 
ated phase fluctuation. The outer scale of turbulence of the model is ~44% of 
A^ATM and at this outer scale of turbulence, given the telescope aperture used 
in the simulation model, the variance of the tip-tilt angle of an image formed 
with the phase fluctuation is reduced by ^40% fr om its expected value base d 
on a KM model which has an infinite outer scale ( Sasiela and Sheltod . 1199,'^ . 



This reduction in tip-tilt fluctuation simulates the active tip-tilt correction at 
SUSI which stabilizes the image of a star. 



7 Testcases 

The main objective of this simulation framework is to test the data reduction 
pipeline of the beam combiners. It is also a good tool to investigate software 
bugs in a data reduction pipeline during its development stage. The following 
sections discuss some testcases using the simulators. 

Several testcases were carried out to demonstrate the functionality of the 
simulator and at the same time the accuracy of the data reduction pipeline for 
both PAVO and MUSCA. Testcase I is carried out to verify the extraction of 
visibility squared, V^, of a set of fringes by the PAVO pipeline. Testcase II- IV 
were carried out to probe the lower bound phase error of the phase-referenced 
fringes constructed by the PAVO and MUSCA phase-referencing pipeline. On 
top of that they are also used to verify the PAVO V"^ pipeline. Table [S] shows 
the input for each testcase. All input parameters for Testcase II-IV are kept 
the same except for the fringe visibility parameter. 
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Fig. 8: Plots of (ZJ^/G.SS)^/^ against r where is either a theoretica l, based 
on a Kolniogorov (KM) model, Eq. (IMl) . or a von Karman (VK) model (| Valley . 
and a simulated structure function. Both axes are normalized to the 



size of the random phases array. The dashed lines in the plots are the differ- 
ence between the solid (KM) and dash dot (simulation) lines relative to the 
former. If the simulated structure function is exactly the same as the theoret- 
ical structure function then the dash lines will be horizontal at 0. The right 
plot is the zoomed in version of the left. 



7.1 To verify the PAVO reduction pipeline 

Testcase TIV are used to verify the reduction pipeline of PAVO at SUSI. 
In these testcases the PAVO simulator is used to investigate the pipeline's 
ability to reproduce the square of the visibility of a given model. The input 
fringe visibility for Testcase I is a sinusoidal function which models a binary 
star system with the primary and secondary stars almost equal in brightness 
(contrast ratio of ^^0.95) and have a projected separation of ~0.04" using a 
15m baseline. On the other hand, the input fringe visibilities for Testcase II- 
IV are constant to model an unresolved single star but have different values 
of instrument visibilities. Since PAVO at SUSI uses only two telescopes at 
any one time, ^ and D each reduces to a 2 x 2 square matrix. The values of 
off-diagonal elements of the matrices, Ci,2 and di,2, are shown in Fig. |9l All 
testcases have 3 sub-cases (A, B and C) where the photon rate is varied to 
simulate an observation of a zeroth, 2nd or 4th magnitude star. The number 
of photons in a generated PAVO frame is adjusted to match the number in 
an actual frame recorded by the PAVO camera at a gain of 5, for a star of 
magnitude mv = 0.0 and my = 2.0, and 25, for a star of magnitude my — 4.0. 

The four types of frames generated by the PAVO simulator are similar to 
those shown in Fig. |31[c) and therefore will not be shown again in this sec- 
tion. The estimated V"^ of each simulated data set extracted by the PAVO 
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Table 5: Input parameters for all testcases 



Testcases — 
Names 4- 



III 



IV 



Section 
^START 

^STEP 

A^S-FITS 
■^R-FITS 
-^Vp-FITS 
^D-FITS 
iVsTEP 

i-SCAN 

A'"med 

c 

-'Vtel 
B 



ray 

V2 

D,d 

ro 
-1 



TO 

Lo 



Section ITU Section 0& EH 
2455848.5 

PAVO: 5.0ms 
MUSCA: 0.2ms 

150 
20 
20 
20 

1024 
150 

- 140/^m 
3 (air, BK7, F7) 
See Fig. |9] 
2 (N3, SI) 
IB1.2I -15m 

Sub-case A; 0.0 
Sub-case B: 2.0 
Sub-case C: 4.0 



0.01-0.36 



0.25 0.10 
See Fig. H 
10cm 
0.5/im 
1ms 
100m 



0.01 




Fig. 9: The value of Ci,2 and c?i_2 used in the simulation versus the frame 
number. The latter is a pseudo-random variable where the difference in value 
between each step is di^2{ta) — ^i, 2(^0-1) ^ A/'(0, O.OI)^m. 
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18 20 22 24 26 28 30 18 20 22 24 26 28 30 



oB, j, (10" rad"') oB, ^ (10* rad"') 

(a) (b) 

Fig. 10: A comparison between user specified models of (solid line) and the 
estimated produced by the PAVO reduction pipeline. The top plot of (a) 
shows the comparison between a wavelength dependent model (input) and 
the (output) values estimated by the pipeline at high photon rate (sub-case A) 
while the bottom plot of (a) represents a similar comparison at lower photon 
rate (sub-case B). The ratio of estimated to model for Testcase I (-I-, x), 
II (o,A), III (V,i>) and IVB (<i) are plotted against spatial frequency in (b). 
Each pair of symbols in the parentheses except for Testcase IV represent sub- 
case A and B respectively. Also in (b), the ratios are shown when atmospheric 
phase noise was excluded in the simulation (no ATM). 



pipeline are plotted against the models fed into the simulator in Fig. [TO] Also 
plotted for comparison, in Fig. [Til are the estimated of Testcase II-IV. 
Both figures show that the estimated values are consistent and within 20% of 
the model. The expected reduction in visibility of the fringes, which is more 
pronounced at higher spatial frequencies (right half of the graphs) , is due to 
the atmospheric phase noise. The observed trend in Fig. ITlTb'). is due to a 
wavelength-dependent scaling factor in the order of unity which is in turn 
caused by the shape of a Fourier domain windowing function used in the 
data reduction pipeline. This scaling factor can be calibrated out because it 
is identical regardless of the visibility function, as seen in Fig. [TUT b) which 
plots the ratios between the estimated and the model V"^ (or TF) across spa- 
tial frequency. Data points with TF significantly larger than 1 is an effect of 
calculating a ratio with the denominator that has a very small value. 

In summary, these testcases show that the reduction pipeline has some 
dependence on the seeing condition, which can be calibrated out, but no mea- 
surable bias as a function of target brightness and the square of the visibility 
of the fringes. 
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oBi 2(10' rad ') 

(a) 




0.00 
0.015 

0.010 

0.005 
0.000 



18 20 22 24 26 28 3018 20 22 24 26 28 30 
0612(10" rad-^) 

(b) 

Fig. 11: Similar to Fig. [TUTa') . the top, middle and bottom plots of (a) and (b) 
show the comparison between a constant V'^ model (input) and the estimated 
(output) of the pipeline for Testcase II, III and IV respectively. The plots on 
the left column are from sub-case A where the photon rate is higher than 
sub-case B which plots are on the right column. The estimated values of 
sub-case B of Testcase IV are very noise but are still unbiased. The plots in 

(a) are simulated with the presence of atmospheric phase noise while plots in 

(b) are simulated without. 
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Fig. 12: The errors of group delay estimates obtained from the PAVO reduction 
pipehne for Testcase II- IV (from top to bottom). The left and right columns 
are errors for sub-cases A and B respectively. The standard deviation of the 
errors, in unit of nanometers, are stated in the plots. The data for the plots 
in (a) and (b) are obtained with and without atmospheric phase noise in the 
simulations. 



7.2 To verify the PAVO and MUSCA phase-referencing pipeline 

In addition to the goal mentioned in previous section, the main aim of Test- 
case II-IV is to test a newly written PAVO and MUSCA phase-referencing 
pipeline and to probe the lower bound phase error of the phase-referenced 
fringes. Both simulators are fed with the same input for this purpose. Some 
inputs which are common to both simulators, e.g. C ^md D for PAVO, are 
resampled at a different rate for MUSCA because the time steps of the two 
simulators are different. Inputs which are not common to both inputs, e.g. are 
defined separately. 

The role of the PAVO data reduction pipeline in these testcases is to provide 
estimates of the group and phase delay of the fringes. The group delay of the 
fringes is defined by the user-specified input variables, D and d. The errors of 
the group delay estimates are shown in Fig. [12] The standard deviation of the 
error of the group delay estimates increases as the visibility of the fringes and 
the photon rate decrease. This is expected as the signal becomes weaker than 
the photon noise and read noise of the camera. Since the phase delay of the 
fringes is estimated from the group delay estimate, residuals of group delay 
errors which are larger than one wavelength of the fringes produce unreliable 
phase delay estimates. Only reliable estimates of phase delay are chosen and 
applied to stabilize the pos ition of a fringe packet generated by the MUSCA 



simulator (jKok et all 120121 ). The percentage of reliable estimates within an 
observation is a function atmospheric seeing, brightness of the target star and 
the visibility of the stellar fringes. 
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scan length (microns) scan lenglh (microns) 



(a) (b) 

Fig. 13: Simulated MUSCA fringes in a waterfall plot format showing the effect 
of (a) without and (b) with phase-referencing. The vertical axis is time and the 
horizontal axis is the optical delay. The top and bottom plots are of different 
scan direction. 



Waterfall plots in Fig. [T3] show the position of the fringe packet through 
time. They clearly show that the fringe packet has a continuously changing 
position and a relatively constant position before and after the phase delay es- 
timates are applied. The phase variation of the phase- referenced fringes across 
multiple scans, as seen in the plot, can be estimate d from the (ICI^) metric. 



which is a measure of the coherence of the fringes (|Kok et ai 120121) , and is 

given as. 



s 



2 



^^^{IC'dJ^Itt^] (36) 



M 



where Nsc is the number of good scans used in the estimation and au is the 
mean wavenumber of MUSCA fringes, which has a typical value of 1.2/j,m~^. 
The factor fj^^'l^ is used to obtain an unbiased estimate of the sample vari- 
ance. Table m shows the value of (|Cp), A^sc and for all testcases. 

For a scan to be considered good, there must be a continuously reliable 
phase delay estimate throughout at least | of the time to make the scan. The 
number of good scans as well as the value of (|Cp) decrease in the same trend 
as the standard deviation of the errors of the group delay estimates increases 
because residuals of the errors which are larger than one wavelength of the 
fringes produce unreliable phase delay estimates. Testcase IVB is an example 
of an extreme case as there are no good scans found. 

The summation of the phase-referenced fringes in all good scans reproduces 
the fringe packet without the distortion introduced by the atmospheric phase 
noise. The effect of coherent summation of the fringes is shown in Fig.[T31 The 
position of the fringe packet within the scan range of the scanning mirror can 
be determined to an uncertainty defined by Sx/\/Nsc- For example, in order to 
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Tabic 6: Results from PAVO and MUSCA phase-referencing algorithm 









with Al M 




w 


ithout ATM 


Testcases J, 


mv — > 


0.0 


2.0 


4.0 


0.0 


2.0 


4.0 






74 


74 


72 


74 


74 


74 


II {V^ = 0.25) 




0.6407 


0.3740 


0.4112 


0.8236 


0.5678 


0.5067 




Sx (nm) 


126 


188 


179 


84 


143 


156 






74 


58 


21 


74 


74 


74 


III (y2 = 0.10) 




0.6387 


0.2630 


0.1352 


0.7582 


0.3752 


0.3708 




Sx (nm) 


127 


220 


275 


100 


188 


189 




Nsc 
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74 








IV (V^ = 0.01) 




0.2684 


N/A 


N/A 


0.3346 


N/A 


N/A 




Sx (nm) 


218 


N/A 


N/A 


198 


N/A 


N/A 




OPD (^m) wavenumbers (1/nm) 



(a) (b) 

Fig. 14: The top and bottom plots of (a) show the result of summation of 
multiple scans of incoherent and coherent fringes respectively. Similarly the 
top and bottom plots of (b) show the huge (about an order of magnitude) 
difference in the power spectrum of the incoherently and coherently summed 
fringe packets respectively. The number of scans used for the summation and 
the (|Cp) metric are indicated in (b). 



determine the fringe packet position to an uncertainty of 5nm, which translates 
to an astrometric uncertainty of 10/xas with 100m baseline, in Testcase IIIB, 
at least 1900 good scans are required. This number is not impra ctical as a 
perio d of good seeing at SUSI can produce over 1000 good scans (|Kok et al . 



20121 ) in about 10 minutes of observation (a total of 3600 scans) with /? Cru, 
which has a V magnitude of 1.3 and an uncalibrated range of 0.02-0.20 as 
measured by PAVO. 

Although the results in Table[n]were simulated for the setup at SUSI, they 
can be extrapolated to estimate the performance of a similar setup at different 
observation sites. For example, the scenario in Testcase IIC (4th magnitude 
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star, seeing condition of ^1" and diameter of light collecting aperture of 14cm) 
is equivalent to observing a ~6th magnitude star at a site with ~0.5" seeing 
(e.g. at NPOil or in Antarctica) but with twice the aperture size (28cm). In 
such scenario, at least 1000 good scans are required to achieve an astrometric 
precision of lO/ias with a 100m baseline. The phase-referencing performances 
in Testcase IIB and IIC are similar because a higher EMCCD gain was used 
to compensate the lower flux in the latter. 

8 Conclusion 

The simulators developed in this work not only have been a very useful tool to 
test the data reduction pipeline developed for the new setup and observation 
technique at SUSI but have also shown the feasibility of visible wavelength 
phase-referencing at sub-wavelength OPD uncertainty. The PAVO-MUSCA 
setup at SUSI can determine the position of a fringe packet of a 2-4th mag- 
nitude star to an accuracy of 5nm by coherently integrating 1000-2000 good 
fringe packet scans. The simulated performance of the dual beam combiner can 
be extrapolated to estimate performance of a similar setup at future possible 
sites (e.g. NPOI and Antarctica). In addition to that, due to the selection of 
input parameters the simulators were designed to accept they could also be 
used to perform simulation for many other functions (e.g. to explore the op- 
tion of expanding the capability of PAVO at CHARA from a 3-telescope to a 
4-telescope beam combiner or to investigate the effect of optical aberration of 
lenses on the performance of the beam combiners) which is beyond its main 
role described in this paper. The IDL code for the simulators can be obtained 
from the corresponding author via email. 
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